Alessandro Filazzola, Batbaatar Amgaa, Charlotte Brown, Issac Heida, Jessica Grenke, Margarete Dettlaff, Tan Bao, & JC Cahill
library(tidyverse)
library(PRISMAstatement)
Conduct a meta-analysis of the literature testing the indirect effects of grazing on animal taxa’s through the direct effects on the plant community.
| date | task |
|---|---|
| Nov 9 | Establish search terms to be used in the meta-analysis |
| Nov 12 | Compile list of journal artcles and sub-divide for each researcher |
| Nov 14 | Begin reviewing papers and extracting data |
| Jan 28 | Complete data extraction from papers |
| Feb 11 | Complete preliminary analysis and set structure for MS |
| Feb 25 | Settle on analyses to be used and begin writing manuscript |
| March 11 | Complete first draft of MS and pass to co-authors |
| March 25 | Comments passed back on draft |
| April 2 | Complete revisions and submit to journal |
A systematic literature search was conducted using Web of Science for all emperical research articles. The review will include all studies globally. The intended purpose of this search is to capture all articles that have documented grazing either along a gradient (e.g. different frequencies or intensity) or presence/absence (e.g. excluded, ungrazed, or retired ranch lands). We condcuted two separate searches to capture studies that tested gradients and studies that compared grazing to ungrazed treatments. Duplicate articles between the searches were removed. We also intentionally excluded terms that resulted in articles not relevant to the purpose of this study including: review, synthesis, policy, social, carbon, and fish. The search terms that used were:
Search A graz* OR livestock AND exclosure* OR exclusion OR exclude* OR ungrazed OR retire* OR fallow*
Search B grazing intensity OR grazing gradient OR stocking rate OR rotation* grazing
This steps includes a. checking for duplicating, b. reviewing each instance for relevancy, c. consistently identifying and documenting exclusion criteria. Outcomes include a list of publications to be used for synthesis, a library of pdfs, and a PRISMA report to ensure the worflow is transparent and reproducible. Papers were excluded with the following characteristics:
evidence <- read.csv("data//ReviewerSearch//evidence.new.csv")
### Identify studies that were excluded
excludes <- evidence %>% group_by(reason.simplified) %>% count(reason.simplified) %>% filter(reason.simplified!="")
## Generate plot
excludes$reason.simplified <- factor(excludes$reason.simplified, levels = excludes$reason.simplified[order(excludes$n)]) ## sort by frequency
ggplot(excludes %>% filter(n > 5) , aes(x=reason.simplified, y=n)) + geom_bar(stat="identity") + coord_flip() + xlab("Reason for exclusion") + theme(text = element_text(size=20))
## frequency of study
year.rate <- evidence %>% group_by(Publication.Year) %>% summarize(n=length(Publication.Year))
ggplot(tail(year.rate,30)) + geom_bar(aes(x=Publication.Year, y=n), stat="identity") + ylab("number of published studies") +xlab("year published") +theme(text = element_text(size=16))
## total number of papers found
nrow(evidence)
## [1] 3488
## number of papers found outside of WoS
other <- read.csv("data/synthesisdata//other.sources.csv")
nrow(other)
## [1] 0
## number of articles excluded
excludes <- evidence %>% filter(exclude=="yes")
nrow(excludes)
## [1] 3149
## relevant papers
review <- evidence %>% filter(exclude!="yes")
nrow(review)
## [1] 339
## papers for meta
datasets <- read.csv("data//binary.simplified.csv")
meta <- length(unique(datasets$uniqueID))
meta
## [1] 239
prisma(found = 3485,
found_other = 4,
no_dupes = 3489,
screened = 3489,
screen_exclusions = 0,
full_text = 3489,
full_text_exclusions = 3250,
qualitative = 239,
quantitative = 90,
width = 800, height = 800)
meta <- read.csv("data//binary.simplified.csv")
## Load packages and functions
library(reshape2)
##
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
##
## smiths
library(metafor)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
## Loading 'metafor' package (version 2.1-0). For an overview
## and introduction to the package please type: help(metafor).
source("meta.evalSE.r") ## Multiple aggregate
## Simplify the estimates column
meta2 <- meta
est <- read.csv("data//Unique.Estimates.Column.csv")
meta2 <- merge(meta2, est, by="Estimate")
## drop Lichen, Microbes and fungi
meta2 <- meta2 %>% filter(Functional.Group != "Producer") %>% filter(Higher.Taxa != "Microbial") ## examine animals only
## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
meta2 <- merge(meta2, ms.data, by="uniqueID")
## subset for abundance only
meta2 <- meta2 %>% filter(est.reduced == "diversity")
meta2[,"functional.taxa"] <- paste0(meta2$Higher.Taxa, meta2$Functional.Group)
## Create Unique identifier column
meta2[,"UniqueSite"] <- paste(meta2$uniqueID, meta2$Higher.Taxa, meta2$Taxa, meta2$functional.taxa, meta2$estimate.simplified, sep="-")
## convert se to sd
meta2[meta2$Stat=="sd","Value"] <- meta2[meta2$Stat=="sd","Value"] / sqrt(meta2[meta2$Stat=="sd","replicate"] )
meta2[meta2$Stat=="sd","Stat"] <- "se"
## drop comparisons that are not pairwise
meta2 <- meta2 %>% filter(grazing.reduced == "ungrazed" | grazing.reduced == "grazed")
meta2 <- meta2 %>% filter(!uniqueID %in% c("30","416")) %>% ## drop study 30 and 416 because anomalous
filter(uniqueID != "111") ## data is missing in study 111
## Drop anomalous studies
meta2 <- meta2 <- meta2 %>% filter(!uniqueID %in% c("654", "1660")) ## native grazer not ungulate
## Test one higher taxa at a time
meta2 <- meta2 %>% filter(Functional.Group != "Other")
## Test domestic grazers only
meta2 <- meta2 %>% filter(grazer.simplified %in% c("domestic","both")) %>% filter(!grazer.spp.reduced == "domestic")
## Use function to extract summary statistics for comparisons
## meta.eval arguments are (meta.data, compare, ids , stats)
grazed.compare <- meta.eval(meta2, grazing.reduced, UniqueSite, Stat)
## Combine the lists into same dataframe
## Rename Columns in second dataframe
grazed.stat <- grazed.compare [[2]] ## extracted statistics
names(grazed.stat) <- c("UniqueSite","grazed_mean","grazed_se","ungrazed_mean","ungrazed_se","grazed_n","ungrazed_n") ## rename columns to match
grazed.raw <- grazed.compare[[1]] ## calculated statistics from raw values
## Join two dataframes
meta.stat <- rbind(grazed.raw, grazed.stat[, names(grazed.raw)])
## Calculate effect size from SE rather than SD. To calculate variance treat replicates as 1.
## https://stats.stackexchange.com/questions/152172/meta-analysis-in-r-with-standard-errors-instead-of-standard-deviations-metafor
meta.stat[,"dummyN"] <- 1
meta.ready <- escalc(n1i = dummyN, n2i = dummyN, m1i = ungrazed_mean, m2i = grazed_mean, sd1i = ungrazed_se, sd2i = grazed_se, data = meta.stat, measure = "ROM", append = TRUE)
## clean up meta.ready
meta.ready <- na.omit(meta.ready) ## drop NAs
## separate out the identifiers
siteID <- matrix(unlist(strsplit(meta.ready$UniqueSite,"-")),ncol=5, byrow=TRUE)
siteID <- data.frame(siteID) ## recreate as dataframe
colnames(siteID) <- c("Study","higher.taxa","taxa","FG","measure") ## add column names
meta.ready <- cbind(data.frame(meta.ready), siteID)
#random-effects meta-analysis for grazed vs ungrazed plots
m1 <- rma(yi=yi, vi=vi, data = meta.ready)
summary(m1)
##
## Random-Effects Model (k = 29; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -4.8711 9.7422 13.7422 16.4066 14.2222
##
## tau^2 (estimated amount of total heterogeneity): 0.0427 (SE = 0.0180)
## tau (square root of estimated tau^2 value): 0.2065
## I^2 (total heterogeneity / total variability): 72.39%
## H^2 (total variability / sampling variability): 3.62
##
## Test for Heterogeneity:
## Q(df = 28) = 109.0878, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1807 0.0504 3.5848 0.0003 0.0819 0.2794 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#mixed-effects meta-analysis for grazed vs ungrazed plots
m2 <- rma(yi=yi, vi=vi, mods=~ -1 + FG, data = subset(meta.ready, vi<4), digits=4 )
summary(m2)
##
## Mixed-Effects Model (k = 29; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -1.2542 2.5085 22.5085 32.4658 46.9529
##
## tau^2 (estimated amount of residual heterogeneity): 0.0248 (SE = 0.0151)
## tau (square root of estimated tau^2 value): 0.1574
## I^2 (residual heterogeneity / unaccounted variability): 56.54%
## H^2 (unaccounted variability / sampling variability): 2.30
##
## Test for Residual Heterogeneity:
## QE(df = 20) = 41.3645, p-val = 0.0033
##
## Test of Moderators (coefficients 1:9):
## QM(df = 9) = 32.6884, p-val = 0.0002
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## FGInvertebrateDetritivores 0.1089 0.0963 1.1306 0.2582 -0.0799 0.2978
## FGInvertebrateHerbivorous 0.3712 0.1242 2.9887 0.0028 0.1278 0.6147
## FGInvertebrateOmnivorous 0.0803 0.2591 0.3098 0.7567 -0.4276 0.5882
## FGInvertebrateParasitic 0.3580 0.2036 1.7584 0.0787 -0.0410 0.7571
## FGInvertebratePollinator 0.3032 0.0772 3.9255 <.0001 0.1518 0.4545
## FGInvertebratePredator 0.2389 0.1384 1.7263 0.0843 -0.0323 0.5101
## FGVertebrateHerbivorous -0.0473 0.1315 -0.3597 0.7190 -0.3050 0.2104
## FGVertebrateOmnivorous -0.0864 0.1175 -0.7349 0.4624 -0.3167 0.1440
## FGVertebratePredator -0.4832 1.0064 -0.4801 0.6312 -2.4557 1.4893
##
## FGInvertebrateDetritivores
## FGInvertebrateHerbivorous **
## FGInvertebrateOmnivorous
## FGInvertebrateParasitic .
## FGInvertebratePollinator ***
## FGInvertebratePredator .
## FGVertebrateHerbivorous
## FGVertebrateOmnivorous
## FGVertebratePredator
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#The benchmark values for I2 are 25, 50 and 75% for low, moderate and high heterogeneity, respectively (Higgins et al. 2003).
## number of studies per comparison
nstudies <- meta.ready %>% group_by(FG) %>% summarise(n=length(FG))
nstudies
## # A tibble: 9 x 2
## FG n
## <fct> <int>
## 1 InvertebrateDetritivores 5
## 2 InvertebrateHerbivorous 4
## 3 InvertebrateOmnivorous 1
## 4 InvertebrateParasitic 1
## 5 InvertebratePollinator 7
## 6 InvertebratePredator 2
## 7 VertebrateHerbivorous 5
## 8 VertebrateOmnivorous 3
## 9 VertebratePredator 1
#write.csv(meta.ready, "data//effectSize//div.animal.csv", row.names=FALSE)
## compare inclusion of moderators
(.9539-.9352)/.9352 ## explains an additional 2%
## [1] 0.01999572
## Produce a forest plot to determine the effect sizes for each study
forest(m2)
regtest(m2)
##
## Regression Test for Funnel Plot Asymmetry
##
## model: mixed-effects meta-regression model
## predictor: standard error
##
## test for funnel plot asymmetry: z = 0.7627, p = 0.4456
confint(m1)
##
## estimate ci.lb ci.ub
## tau^2 0.0427 0.0149 0.0800
## tau 0.2065 0.1220 0.2829
## I^2(%) 72.3950 47.7645 83.1039
## H^2 3.6225 1.9144 5.9185
## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m2, ylim=c(1.2,0))
## Calculate rosenthals Failsafe number
fsn(yi, vi, data=meta.ready)
##
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 384
# ## generate plot with spaces inbetween
# forest(m1, atransf=exp, cex=0.75, ylim=c(-1, 24),
# order=order(meta.ready$GI.compare,meta.ready$taxa), rows=c(3:4,7,10:16,19:21),
# # mlab="RE model for all studies", psize=1, slab= paste(meta.ready$Study, meta.ready$taxa, meta.ready$measure))
#
# addpoly(res.w, row=18, cex=0.75, atransf=exp, mlab="RE model for green wall")
# addpoly(res.r, row= 9, cex=0.75, atransf=exp, mlab="RE model for green roof")
# addpoly(res.rd, row= 6, cex=0.75, atransf=exp, mlab="RE model for roadsides")
# addpoly(res.p, row= 2, cex=0.75, atransf=exp, mlab="RE model for retention ponds")
meta <- read.csv("data//binary.simplified.csv")
## Load packages and functions
library(reshape2)
library(metafor)
source("meta.evalSE.r") ## Multiple aggregate
## Simplify the estimates column
meta2 <- meta
est <- read.csv("data//Unique.Estimates.Column.csv")
meta2 <- merge(meta2, est, by="Estimate")
## drop duplicates
ests <- meta2[!duplicated(meta2[,c(1,21:22)]),c(1,21:22)]
ests.used <- ests %>% filter(!est.reduced %in% c("omit", "behavioural","fitness")) ## filter only the variables used
## write estimates
write.csv(ests.used, "data//appendixA.csv", row.names = F)
meta.ready <- read.csv("abundanceEff.csv")
#random-effects meta-analysis for grazed vs ungrazed plots
m1 <- rma(yi=yi, vi=vi, data = meta.ready, test="t" )
summary(m1)
##
## Random-Effects Model (k = 122; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -140.5963 281.1926 285.1926 290.7842 285.2943
##
## tau^2 (estimated amount of total heterogeneity): 0.4202 (SE = 0.0719)
## tau (square root of estimated tau^2 value): 0.6482
## I^2 (total heterogeneity / total variability): 92.78%
## H^2 (total variability / sampling variability): 13.85
##
## Test for Heterogeneity:
## Q(df = 121) = 2302.4568, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.2089 0.0692 3.0203 0.0031 0.0720 0.3459 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#mixed-effects meta-analysis for grazed vs ungrazed plots
m2 <- rma(yi=yi, vi=vi, mods=~ -1 + FG, data = subset(meta.ready, vi<4), digits=4, test="t" )
summary(m2)
##
## Mixed-Effects Model (k = 122; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -127.2609 254.5219 274.5219 301.7958 276.6788
##
## tau^2 (estimated amount of residual heterogeneity): 0.3838 (SE = 0.0696)
## tau (square root of estimated tau^2 value): 0.6195
## I^2 (residual heterogeneity / unaccounted variability): 89.86%
## H^2 (unaccounted variability / sampling variability): 9.86
##
## Test for Residual Heterogeneity:
## QE(df = 113) = 1723.5800, p-val < .0001
##
## Test of Moderators (coefficients 1:9):
## F(df1 = 9, df2 = 113) = 3.0013, p-val = 0.0030
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## FGInvertebrateDetritivores -0.3829 0.1930 -1.9839 0.0497 -0.7654 -0.0005
## FGInvertebrateHerbivorous -0.0423 0.1854 -0.2284 0.8198 -0.4097 0.3250
## FGInvertebrateOmnivorous 0.3817 0.3808 1.0025 0.3183 -0.3727 1.1362
## FGInvertebrateParasitic 0.1809 0.2729 0.6630 0.5087 -0.3597 0.7215
## FGInvertebratePollinator 0.2097 0.2163 0.9697 0.3343 -0.2188 0.6382
## FGInvertebratePredator 0.1892 0.2211 0.8557 0.3940 -0.2489 0.6273
## FGVertebrateHerbivorous 0.4626 0.1370 3.3771 0.0010 0.1912 0.7340
## FGVertebrateOmnivorous 0.2931 0.1947 1.5055 0.1350 -0.0926 0.6788
## FGVertebratePredator 0.4841 0.1939 2.4971 0.0140 0.1000 0.8682
##
## FGInvertebrateDetritivores *
## FGInvertebrateHerbivorous
## FGInvertebrateOmnivorous
## FGInvertebrateParasitic
## FGInvertebratePollinator
## FGInvertebratePredator
## FGVertebrateHerbivorous **
## FGVertebrateOmnivorous
## FGVertebratePredator *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(snow)
hetero <- gosh(m1, subsets=10000, parallel = "snow", ncpus=40)
## Fitting 10000 models (based on random subsets).
plot(hetero, breaks=100)
### Bias tests
regtest(m2)
##
## Regression Test for Funnel Plot Asymmetry
##
## model: mixed-effects meta-regression model
## predictor: standard error
##
## test for funnel plot asymmetry: t = -1.4297, df = 112, p = 0.1556
## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m2, ylim=c(1.2,0))
## Calculate rosenthals Failsafe number
fsn(yi, vi, data=meta.ready)
##
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 8156
meta.ready <- read.csv("diversityEff.csv")
#random-effects meta-analysis for grazed vs ungrazed plots
m1 <- rma(yi=yi, vi=vi, data = meta.ready)
summary(m1)
##
## Random-Effects Model (k = 29; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -4.8711 9.7422 13.7422 16.4066 14.2222
##
## tau^2 (estimated amount of total heterogeneity): 0.0427 (SE = 0.0180)
## tau (square root of estimated tau^2 value): 0.2065
## I^2 (total heterogeneity / total variability): 72.39%
## H^2 (total variability / sampling variability): 3.62
##
## Test for Heterogeneity:
## Q(df = 28) = 109.0878, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1807 0.0504 3.5848 0.0003 0.0819 0.2794 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hetero <- gosh(m1, subsets=10000, parallel = "snow", ncpus=40)
## Fitting 10000 models (based on random subsets).
plot(hetero, breaks=100)
#mixed-effects meta-analysis for grazed vs ungrazed plots
m2 <- rma(yi=yi, vi=vi, mods=~ -1 + FG, data = subset(meta.ready, vi<4), digits=4, test="t" )
summary(m2)
##
## Mixed-Effects Model (k = 29; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -1.2542 2.5085 22.5085 32.4658 46.9529
##
## tau^2 (estimated amount of residual heterogeneity): 0.0248 (SE = 0.0151)
## tau (square root of estimated tau^2 value): 0.1574
## I^2 (residual heterogeneity / unaccounted variability): 56.54%
## H^2 (unaccounted variability / sampling variability): 2.30
##
## Test for Residual Heterogeneity:
## QE(df = 20) = 41.3645, p-val = 0.0033
##
## Test of Moderators (coefficients 1:9):
## F(df1 = 9, df2 = 20) = 3.6320, p-val = 0.0078
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## FGInvertebrateDetritivores 0.1089 0.0963 1.1306 0.2716 -0.0920 0.3099
## FGInvertebrateHerbivorous 0.3712 0.1242 2.9887 0.0073 0.1121 0.6303
## FGInvertebrateOmnivorous 0.0803 0.2591 0.3098 0.7599 -0.4603 0.6208
## FGInvertebrateParasitic 0.3580 0.2036 1.7584 0.0940 -0.0667 0.7828
## FGInvertebratePollinator 0.3032 0.0772 3.9255 0.0008 0.1421 0.4643
## FGInvertebratePredator 0.2389 0.1384 1.7263 0.0997 -0.0498 0.5276
## FGVertebrateHerbivorous -0.0473 0.1315 -0.3597 0.7228 -0.3216 0.2270
## FGVertebrateOmnivorous -0.0864 0.1175 -0.7349 0.4709 -0.3316 0.1588
## FGVertebratePredator -0.4832 1.0064 -0.4801 0.6364 -2.5825 1.6161
##
## FGInvertebrateDetritivores
## FGInvertebrateHerbivorous **
## FGInvertebrateOmnivorous
## FGInvertebrateParasitic .
## FGInvertebratePollinator ***
## FGInvertebratePredator .
## FGVertebrateHerbivorous
## FGVertebrateOmnivorous
## FGVertebratePredator
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Bias tests
regtest(m2)
##
## Regression Test for Funnel Plot Asymmetry
##
## model: mixed-effects meta-regression model
## predictor: standard error
##
## test for funnel plot asymmetry: t = 0.7627, df = 19, p = 0.4550
## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m2, ylim=c(0.8,0))
## Calculate rosenthals Failsafe number
fsn(yi, vi, data=meta.ready)
##
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 384
meta <- read.csv("data//plantAnalysis.csv")
meta2 <- meta
## Load packages and functions
library(reshape2)
library(metafor)
source("meta.evalSE.r") ## Multiple aggregate
## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
meta2 <- merge(meta2, ms.data, by="uniqueID")
## Drop mychoriza because only two studies
## Create Unique identifier column
meta2 <- meta2 %>% filter(estimate.simplified != "omit")
meta2[,"UniqueSite"] <- paste(meta2$uniqueID, meta2$Taxa, meta2$estimate.simplified, sep="-")
## convert se to sd
meta2[meta2$Stat=="sd","Value"] <- meta2[meta2$Stat=="sd","Value"] / sqrt(meta2[meta2$Stat=="sd","replicate"] )
meta2[meta2$Stat=="sd","Stat"] <- "se"
## drop comparisons that are not pairwise
meta2 <- meta2 %>% filter(grazing.reduced == "ungrazed" | grazing.reduced == "grazed")
meta2 <- meta2 %>% filter(uniqueID != "30") ## drop study 30 because anomalous
## Drop anomalous studies
meta2 <- meta2 <- meta2 %>% filter(!uniqueID %in% c("654", "1660")) ## native grazer not ungulate
meta2 <- meta2 %>% filter(estimate.simplified != "fitness") ## drop fitness because only study
## Test domestic grazers only
meta2 <- meta2 %>% filter(grazer.simplified %in% c("domestic","both")) %>%
filter(grazer.spp.reduced != "domestic") ## select only sheep and cattle comparisons
## Use function to extract summary statistics for comparisons
## meta.eval arguments are (meta.data, compare, ids , stats)
grazed.compare <- meta.eval(meta2, grazing.reduced, UniqueSite, Stat)
## Combine the lists into same dataframe
## Rename Columns in second dataframe
grazed.stat <- grazed.compare [[2]] ## extracted statistics
names(grazed.stat) <- c("UniqueSite","grazed_mean","grazed_se","ungrazed_mean","ungrazed_se","grazed_n","ungrazed_n") ## rename columns to match
grazed.raw <- grazed.compare[[1]] ## calculated statistics from raw values
## Join two dataframes
meta.stat <- rbind(grazed.raw, grazed.stat[, names(grazed.raw)])
## Calculate effect size from SE rather than SD. To calculate variance treat replicates as 1.
## https://stats.stackexchange.com/questions/152172/meta-analysis-in-r-with-standard-errors-instead-of-standard-deviations-metafor
meta.stat[,"dummyN"] <- 1
meta.ready <- escalc(n1i = dummyN, n2i = dummyN, m1i = ungrazed_mean, m2i = grazed_mean, sd1i = ungrazed_se, sd2i = grazed_se, data = meta.stat, measure = "ROM", append = TRUE)
## Warning in escalc(n1i = dummyN, n2i = dummyN, m1i = ungrazed_mean, m2i =
## grazed_mean, : Some 'yi' and/or 'vi' values equal to +-Inf. Recoded to NAs.
## clean up meta.ready
meta.ready <- na.omit(meta.ready) ## drop NAs
## separate out the identifiers
siteID <- matrix(unlist(strsplit(meta.ready$UniqueSite,"-")),ncol=3, byrow=TRUE)
siteID <- data.frame(siteID) ## recreate as dataframe
colnames(siteID) <- c("Study","taxa","measure") ## add column names
meta.ready <- cbind(data.frame(meta.ready), siteID)
#random-effects meta-analysis for grazed vs ungrazed plots
m1 <- rma(yi=yi, vi=vi, data = meta.ready)
## Warning in rma(yi = yi, vi = vi, data = meta.ready): There are outcomes with
## non-positive sampling variances.
## Warning in rma(yi = yi, vi = vi, data = meta.ready): Cannot compute Q-test, I^2,
## or H^2 with non-positive sampling variances.
summary(m1)
##
## Random-Effects Model (k = 65; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -60.9895 121.9791 125.9791 130.2968 126.1758
##
## tau^2 (estimated amount of total heterogeneity): 0.3054 (SE = 0.0621)
## tau (square root of estimated tau^2 value): 0.5527
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2942 0.0741 3.9692 <.0001 0.1489 0.4395 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## get last year grazed as covariate
yr.grazed <- meta2 %>% group_by(uniqueID) %>% summarize(yr=max(yr.grazed))
colnames(yr.grazed)[1] <- "Study"
site.yr <- merge(siteID, yr.grazed, by="Study")
meta.ready[,"yr.grazed"] <- site.yr$yr
meta.ready[,"previously.grazed"] <- as.factor(ifelse(is.na(meta.ready$yr.grazed),0,1))
#mixed-effects meta-analysis for grazed vs ungrazed plots
m2 <- rma(yi=yi, vi=vi, mods=~ -1 + measure , data = subset(meta.ready, vi!=0))
summary(m2)
##
## Mixed-Effects Model (k = 61; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -50.3340 100.6679 110.6679 120.8832 111.8444
##
## tau^2 (estimated amount of residual heterogeneity): 0.2471 (SE = 0.0553)
## tau (square root of estimated tau^2 value): 0.4970
## I^2 (residual heterogeneity / unaccounted variability): 95.65%
## H^2 (unaccounted variability / sampling variability): 23.00
##
## Test for Residual Heterogeneity:
## QE(df = 57) = 485.3761, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## QM(df = 4) = 36.9574, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## measureabundance 0.4267 0.0872 4.8932 <.0001 0.2558 0.5976
## measuredebris 0.9367 0.2710 3.4568 0.0005 0.4056 1.4678
## measurediversity -0.1251 0.1415 -0.8839 0.3768 -0.4024 0.1523
## measureMycorrhizaColonize 0.1887 0.3545 0.5324 0.5944 -0.5060 0.8835
##
## measureabundance ***
## measuredebris ***
## measurediversity
## measureMycorrhizaColonize
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## number of studies
meta.ready %>% group_by(measure) %>% summarize(n=length(measure))
## # A tibble: 4 x 2
## measure n
## <fct> <int>
## 1 abundance 41
## 2 debris 5
## 3 diversity 17
## 4 MycorrhizaColonize 2
fsn(yi, vi, data=meta.ready)
##
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: NA
## Target Significance Level: 0.05
##
## Fail-safe N: NaN
## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m2)
regtest(m2)
##
## Regression Test for Funnel Plot Asymmetry
##
## model: mixed-effects meta-regression model
## predictor: standard error
##
## test for funnel plot asymmetry: z = 0.9919, p = 0.3212
## write
write.csv(subset(meta.ready, vi<4), "plantEff.csv", row.names=F)
meta.ready <- read.csv("plantEff.csv")
## Abundance
#random-effects meta-analysis for grazed vs ungrazed plots
m1 <- rma(yi=yi, vi=vi, data = subset(meta.ready, measure == "abundance" & vi!=0))
summary(m1)
##
## Random-Effects Model (k = 40; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -32.6182 65.2364 69.2364 72.5635 69.5697
##
## tau^2 (estimated amount of total heterogeneity): 0.2399 (SE = 0.0659)
## tau (square root of estimated tau^2 value): 0.4898
## I^2 (total heterogeneity / total variability): 96.06%
## H^2 (total variability / sampling variability): 25.38
##
## Test for Heterogeneity:
## Q(df = 39) = 340.6311, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.4261 0.0861 4.9478 <.0001 0.2573 0.5949 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hetero <- gosh(m1, subsets=10000, parallel = "snow", ncpus=10)
## Fitting 10000 models (based on random subsets).
plot(hetero, breaks=100)
## Diversity
#random-effects meta-analysis for grazed vs ungrazed plots
m2 <- rma(yi=yi, vi=vi, data = subset(meta.ready, measure == "diversity" & vi!=0))
summary(m2)
##
## Random-Effects Model (k = 14; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -0.9716 1.9432 5.9432 7.0731 7.1432
##
## tau^2 (estimated amount of total heterogeneity): 0.0321 (SE = 0.0208)
## tau (square root of estimated tau^2 value): 0.1793
## I^2 (total heterogeneity / total variability): 67.07%
## H^2 (total variability / sampling variability): 3.04
##
## Test for Heterogeneity:
## Q(df = 13) = 34.4838, p-val = 0.0010
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.1123 0.0630 -1.7805 0.0750 -0.2358 0.0113 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
hetero <- gosh(m2, subsets=10000, parallel = "snow", ncpus=10)
## Fitting 10000 models (based on random subsets).
plot(hetero, breaks=100)
## load metadata of manuscripts
ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
ms.data[is.na(ms.data)] <- 0
## Warning in `[<-.factor`(`*tmp*`, thisvar, value = 0): invalid factor level, NA
## generated
## Warning in `[<-.factor`(`*tmp*`, thisvar, value = 0): invalid factor level, NA
## generated
## load effect sizes from each model
div <- read.csv("diversityEff.csv")
abd <- read.csv("abundanceEff.csv")
plt <- read.csv("plantEff.csv")
## list unique studies
studies <- c(div$Study, abd$Study, plt$Study)
studies <- studies[!duplicated(studies)]
## select only used studies
msUsed <- ms.data[ms.data$uniqueID %in% studies,]
## Which grazer
msUsed %>% group_by(grazer.spp.reduced) %>% summarize(n=length(grazer.spp.reduced))
## # A tibble: 6 x 2
## grazer.spp.reduced n
## <fct> <int>
## 1 both 17
## 2 cattle 58
## 3 cattle and domestic 5
## 4 cattle and indigenous 7
## 5 sheep 20
## 6 sheep and domestic 2
##indigenous grazers
msUsed %>% group_by(indigenous.simplified
) %>% summarize(n=length(indigenous.simplified))
## # A tibble: 8 x 2
## indigenous.simplified n
## <fct> <int>
## 1 0 1
## 2 co-occurrence 1
## 3 co-occurring 33
## 4 historic 15
## 5 none 42
## 6 UNKN 1
## 7 unknown 15
## 8 yes 1
## Compare grazer species on effect size
## Plant abundance and diversity
names(plt)[11] <- "uniqueID"
effdata <- merge(plt, msUsed, by="uniqueID")
m1 <- rma(yi=yi, vi=vi, mods=~ grazer.spp.reduced , data = subset(effdata, vi!=0 & measure == "abundance"))
## Warning in rma(yi = yi, vi = vi, mods = ~grazer.spp.reduced, data =
## subset(effdata, : Redundant predictors dropped from the model.
summary(m1)
##
## Mixed-Effects Model (k = 40; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -29.4268 58.8535 70.8535 80.1856 73.8535
##
## tau^2 (estimated amount of residual heterogeneity): 0.2391 (SE = 0.0698)
## tau (square root of estimated tau^2 value): 0.4889
## I^2 (residual heterogeneity / unaccounted variability): 94.08%
## H^2 (unaccounted variability / sampling variability): 16.90
## R^2 (amount of heterogeneity accounted for): 0.34%
##
## Test for Residual Heterogeneity:
## QE(df = 35) = 292.8272, p-val < .0001
##
## Test of Moderators (coefficients 2:5):
## QM(df = 4) = 3.9913, p-val = 0.4072
##
## Model Results:
##
## estimate se zval pval
## intrcpt 0.7688 0.2249 3.4188 0.0006
## grazer.spp.reducedcattle -0.4362 0.2530 -1.7242 0.0847
## grazer.spp.reducedcattle and domestic -0.2374 0.3605 -0.6585 0.5102
## grazer.spp.reducedcattle and indigenous -0.6543 0.4173 -1.5682 0.1168
## grazer.spp.reducedsheep -0.2753 0.3157 -0.8720 0.3832
## ci.lb ci.ub
## intrcpt 0.3280 1.2095 ***
## grazer.spp.reducedcattle -0.9321 0.0597 .
## grazer.spp.reducedcattle and domestic -0.9441 0.4692
## grazer.spp.reducedcattle and indigenous -1.4721 0.1635
## grazer.spp.reducedsheep -0.8941 0.3435
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m2 <- rma(yi=yi, vi=vi, mods=~ grazer.spp.reduced , data = subset(effdata, vi!=0 & measure == "diversity"))
## Warning in rma(yi = yi, vi = vi, mods = ~grazer.spp.reduced, data =
## subset(effdata, : Redundant predictors dropped from the model.
summary(m2)
##
## Mixed-Effects Model (k = 14; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -1.2390 2.4780 14.4780 15.6613 56.4780
##
## tau^2 (estimated amount of residual heterogeneity): 0.0455 (SE = 0.0342)
## tau (square root of estimated tau^2 value): 0.2133
## I^2 (residual heterogeneity / unaccounted variability): 65.69%
## H^2 (unaccounted variability / sampling variability): 2.91
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 9) = 26.2830, p-val = 0.0018
##
## Test of Moderators (coefficients 2:5):
## QM(df = 4) = 2.9343, p-val = 0.5689
##
## Model Results:
##
## estimate se zval pval
## intrcpt -0.1679 0.1125 -1.4925 0.1356
## grazer.spp.reducedcattle -0.0451 0.1776 -0.2537 0.7997
## grazer.spp.reducedcattle and domestic 0.0694 0.2580 0.2689 0.7880
## grazer.spp.reducedcattle and indigenous 0.3359 0.2260 1.4865 0.1371
## grazer.spp.reducedsheep 0.1247 0.2402 0.5193 0.6036
## ci.lb ci.ub
## intrcpt -0.3884 0.0526
## grazer.spp.reducedcattle -0.3932 0.3030
## grazer.spp.reducedcattle and domestic -0.4362 0.5749
## grazer.spp.reducedcattle and indigenous -0.1070 0.7787
## grazer.spp.reducedsheep -0.3460 0.5954
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Animal abundance
names(abd)[11] <- "uniqueID"
effdata <- merge(abd, msUsed, by="uniqueID")
m1 <- rma(yi=yi, vi=vi, mods=~ grazer.spp.reduced , data = subset(effdata, vi!=0))
## Warning in rma(yi = yi, vi = vi, mods = ~grazer.spp.reduced, data =
## subset(effdata, : Redundant predictors dropped from the model.
summary(m1)
##
## Mixed-Effects Model (k = 122; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -135.2740 270.5480 284.5480 303.8231 285.5850
##
## tau^2 (estimated amount of residual heterogeneity): 0.4121 (SE = 0.0728)
## tau (square root of estimated tau^2 value): 0.6420
## I^2 (residual heterogeneity / unaccounted variability): 91.70%
## H^2 (unaccounted variability / sampling variability): 12.05
## R^2 (amount of heterogeneity accounted for): 1.92%
##
## Test for Residual Heterogeneity:
## QE(df = 116) = 1529.5893, p-val < .0001
##
## Test of Moderators (coefficients 2:6):
## QM(df = 5) = 4.9338, p-val = 0.4240
##
## Model Results:
##
## estimate se zval pval
## intrcpt 0.1559 0.2081 0.7491 0.4538
## grazer.spp.reducedcattle 0.1638 0.2289 0.7157 0.4742
## grazer.spp.reducedcattle and domestic -0.0439 0.5136 -0.0855 0.9318
## grazer.spp.reducedcattle and indigenous 0.2268 0.3348 0.6776 0.4981
## grazer.spp.reducedsheep -0.1416 0.2513 -0.5633 0.5733
## grazer.spp.reducedsheep and domestic -0.2950 0.3885 -0.7593 0.4477
## ci.lb ci.ub
## intrcpt -0.2520 0.5638
## grazer.spp.reducedcattle -0.2848 0.6123
## grazer.spp.reducedcattle and domestic -1.0505 0.9627
## grazer.spp.reducedcattle and indigenous -0.4293 0.8830
## grazer.spp.reducedsheep -0.6342 0.3510
## grazer.spp.reducedsheep and domestic -1.0564 0.4665
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Animal diversity
names(div)[11] <- "uniqueID"
effdata <- merge(div, msUsed, by="uniqueID")
m2 <- rma(yi=yi, vi=vi, mods=~ grazer.spp.reduced , data = subset(effdata, vi!=0 & measure == "diversity"))
## Warning in rma(yi = yi, vi = vi, mods = ~grazer.spp.reduced, data =
## subset(effdata, : Redundant predictors dropped from the model.
summary(m2)
##
## Mixed-Effects Model (k = 3; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -1.0752 2.1504 8.1504 2.1504 32.1504
##
## tau^2 (estimated amount of residual heterogeneity): 0 (SE = 1.7193)
## tau (square root of estimated tau^2 value): 0
## I^2 (residual heterogeneity / unaccounted variability): 0.00%
## H^2 (unaccounted variability / sampling variability): 1.00
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 1) = 0.1172, p-val = 0.7321
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.3649, p-val = 0.5458
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.2618 0.1256 -2.0841 0.0372 -0.5081 -0.0156
## grazer.spp.reducedcattle 0.1469 0.2432 0.6041 0.5458 -0.3297 0.6236
##
## intrcpt *
## grazer.spp.reducedcattle
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
require(ggmap)
require(maps)
### Start with base map of world
mp <- NULL
mapWorld <- borders("world", colour="gray50", fill="gray65") # create a layer of borders
mp <- ggplot() + theme_bw()+ mapWorld
## colorblind-friendly palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#000000")
msUsed <- msUsed %>% filter(lat != "") ## remove blanks for lat and lon
msUsed$lat <- as.numeric(levels(msUsed$lat))[msUsed$lat] ## convert from factor to number
msUsed$lon <- as.numeric(levels(msUsed$lon))[msUsed$lon] ## convert from factor to number
## plot points on top
mp <- mp+ geom_point(data=msUsed , aes(x=lon, y=lat), size=2)
mp
unique(msUsed$grazer.spp.reduced)
## [1] cattle sheep both
## [4] cattle and domestic sheep and domestic cattle and indigenous
## 9 Levels: both cattle cattle and domestic cattle and indigenous ... sheep and indigenous
msUsed[msUsed$grazer.spp.reduced=="domestic",]
## [1] uniqueID Publication.Year StudyYear
## [4] grazer.spp grazer.spp.reduced grazer.status
## [7] grazer.simplified contrast last.grazing.event
## [10] study.duration indigenous.grazers indigenous.simplified
## [13] grazing.frequency lat lon
## [16] elevation country n.sites
## [19] survey.technique
## <0 rows> (or 0-length row.names)
## figure out country from GPS
names <- map.where(database="world", msUsed$lon, msUsed$lat)
## separate subcountry identifiers
country <- gsub("?:.*", "", names)
length(unique(country))
## [1] 27
table(country)
## country
## Argentina Australia Canada Chile China Denmark
## 4 10 5 1 4 1
## France Germany Israel Italy Kazakhstan Kenya
## 1 5 1 1 1 9
## Lesotho Mongolia Netherlands New Zealand Norway Portugal
## 1 1 3 1 1 1
## Romania Senegal South Africa Spain Sweden Switzerland
## 1 1 1 1 3 1
## UK USA
## 13 24
## extract climate
koppen <- read.table("data//KoppenClimate.txt", header=T)
koppen[,"location"] <- paste(koppen$Lat, koppen$Lon, sep="-")
msUsed[,"location"] <- paste(msUsed$lat, msUsed$lon, sep="-")
ClosestMatch2 = function(string, stringVector){
stringVector[stringdist::amatch(string, stringVector, maxDist=Inf)]
}
climID <- ClosestMatch2(msUsed$location, koppen$location)
climNames <- merge(data.frame(location=climID), koppen, by.x="location")
## list climates
climNames %>% group_by(Cls) %>% summarize(n=length(Cls))
## # A tibble: 17 x 2
## Cls n
## <fct> <int>
## 1 Am 2
## 2 Aw 2
## 3 BSh 2
## 4 BSk 20
## 5 BWh 1
## 6 BWk 1
## 7 Cfa 8
## 8 Cfb 24
## 9 Cfc 1
## 10 Csa 6
## 11 Csb 7
## 12 Cwb 2
## 13 Dfb 10
## 14 Dfc 10
## 15 Dwa 3
## 16 Dwb 1
## 17 EF 1
climNames[,"location"] <- as.character(climNames[,"location"] )
names(climNames)[2:3] <- c("lat","lon")
## Drop anomalous point
climNames <- subset(climNames, lat > -80)
## More broad groupings
climNames["group"] <- substr(climNames$Cls, 1, 1)
climgroups <- data.frame(group=c("A","B","C","D"), climate=c("Tropical","Arid/Semi-arid","Temperate","Continental"))
climNames <- merge(climNames, climgroups, by="group")
cbPalette <- c( "#D55E00", "#56B4E9", "#0072B2", "#009E73", "#F0E442", "#E69F00", "#CC79A7","#000000","#999999")
## map with climate
mp + geom_point(data=msUsed , aes(x=lon, y=lat, col=grazer.spp.reduced ), size=2) + xlab("Longitude") + ylab("Latitude") + xlim(-180, 180) + ylim(-90,90)+ scale_color_manual(values = cbPalette, name = "Grazer composition")
## Load in effect sizes
occurdat<-list.files( pattern="Eff.csv$",full=T)
animal.abd <- read.csv(occurdat[2])
# animal.abd <- subset(animal.abd, measure=="diversity" & vi!=0) ## one measure at a time for plants
colnames(animal.abd)[colnames(animal.abd)=="Study"] <- "uniqueID"
## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
meta3 <- merge(animal.abd, ms.data, by="uniqueID")
meta3[,c("yi")] <- meta3[,c("yi")] ## inverse effect
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:tidyr':
##
## extract
#
# ## Convert lat lon to spatial points
# meta3 <- meta3 %>% filter(lat != "") ## remove blanks for lat and lon
# meta3 <- meta3[!grepl(";", meta3$lat),] ## remove multiple entries
# meta3$lat <- as.numeric(levels(meta3$lat))[meta3$lat] ## convert from factor to number
# meta3$lon <- as.numeric(levels(meta3$lon))[meta3$lon] ## convert from factor to number
#
# ## generate dataframe
# gps <- data.frame(uniqueID = meta3$uniqueID, lon=meta3$lon, lat=meta3$lat)
# gps <- gps[!duplicated(gps$uniqueID),]
#
#
# ## assign spatial points
# crs.world <-CRS("+proj=longlat +datum=WGS84")
# coordinates(gps) <- ~lon+lat
# proj4string(gps) <- crs.world
#
# ## load rasters
# temp <- raster("C:\\Users\\Fitz\\Downloads\\bioclim2\\bio01.tif")
# prec <- raster("C:\\Users\\Fitz\\Downloads\\bioclim2\\bio12.tif")
#
# temp2 <- temp
#
# clim <- stack(temp,prec)
#
#
# ## extract temperature and preciptiation
# gps.clim <- data.frame(raster::extract(clim, gps))
# gps.clim["uniqueID"] <- gps$uniqueID
## get CRU data
CRU <- read.csv("data//climateCRU.csv")
CRU <- na.omit(CRU)
## compare climate & grazing duration variables against effect size
## get last year grazed as covariate
meta.ready <- merge(meta3, CRU, by="uniqueID")
### Temperature only
#mixed-effects meta-analysis for grazed vs ungrazed plots
m3 <- rma(yi=yi, vi=vi, mods=~ poly(temp,2), data = meta.ready , test="t" )
summary(m3)
##
## Mixed-Effects Model (k = 27; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -4.1041 8.2082 16.2082 20.9204 18.3135
##
## tau^2 (estimated amount of residual heterogeneity): 0.0377 (SE = 0.0182)
## tau (square root of estimated tau^2 value): 0.1942
## I^2 (residual heterogeneity / unaccounted variability): 68.63%
## H^2 (unaccounted variability / sampling variability): 3.19
## R^2 (amount of heterogeneity accounted for): 16.62%
##
## Test for Residual Heterogeneity:
## QE(df = 24) = 68.9193, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 24) = 2.2044, p-val = 0.1322
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.1789 0.0503 3.5571 0.0016 0.0751 0.2827 **
## poly(temp, 2)1 -0.2404 0.2547 -0.9435 0.3548 -0.7661 0.2854
## poly(temp, 2)2 -0.4595 0.2378 -1.9320 0.0653 -0.9504 0.0314 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## line
preds <- predict(m3, temp=meta.ready$temp)
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#D55E00", "#0072B2","#000000")
## generate column for FG
meta.ready[,"trophic"] <- gsub("^.*?ertebrate","",meta.ready$FG)
### Plot against temperature
ggplot(meta.ready, aes(x=temp, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2)+ geom_jitter(aes(color=trophic, shape=higher.taxa), size=4) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Mean annual temperature (°C)") + ylab("Effect of livestock exclusion (LRR)") + geom_line(data=data.frame(preds), aes(x=meta.ready$temp, y=pred), lwd=2) +
geom_ribbon(data=data.frame(preds), aes(x=meta.ready$temp, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) + scale_color_manual(values=cbPalette)+xlim(0,20) + annotate("text", x=2, y=-.6, label="animal \n diversity", size=6)
## Warning: Ignoring unknown aesthetics: y
### Preciptiation
test <- subset(meta.ready, FG == "VertebrateHerbivorous" | FG == "InvertebrateHerbivorous")
#mixed-effects meta-analysis for grazed vs ungrazed plots
m4 <- rma(yi=yi, vi=vi, mods=~ poly(precip,2), data = subset(meta.ready))
summary(m4)
##
## Mixed-Effects Model (k = 27; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -5.2128 10.4255 18.4255 23.1377 20.5308
##
## tau^2 (estimated amount of residual heterogeneity): 0.0448 (SE = 0.0203)
## tau (square root of estimated tau^2 value): 0.2116
## I^2 (residual heterogeneity / unaccounted variability): 73.22%
## H^2 (unaccounted variability / sampling variability): 3.73
## R^2 (amount of heterogeneity accounted for): 1.05%
##
## Test for Residual Heterogeneity:
## QE(df = 24) = 86.7515, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 1.4233, p-val = 0.4908
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.1685 0.0537 3.1358 0.0017 0.0632 0.2738 **
## poly(precip, 2)1 0.1461 0.2857 0.5112 0.6092 -0.4139 0.7060
## poly(precip, 2)2 -0.2959 0.2836 -1.0433 0.2968 -0.8517 0.2600
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## line
preds <- predict(m4, precip=meta.ready$precip)
### Plot against temperature
ggplot(meta.ready, aes(x=precip, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2)+ geom_jitter(aes(color=trophic, shape=higher.taxa),size=4) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Mean annual preciptiation (mm)") + ylab("Effect of livestock exclusion (LRR)") + scale_color_manual(values=cbPalette) +
geom_line(data=data.frame(preds), aes(x=meta.ready$precip, y=pred), lwd=2) +
geom_ribbon(data=data.frame(preds), aes(x=meta.ready$precip, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2)
## Warning: Ignoring unknown aesthetics: y
## Time since grazed Grazed
## Load in effect sizes
occurdat<-list.files( pattern="Eff.csv$",full=T)
animal.abd <- read.csv(occurdat[1])
colnames(animal.abd)[colnames(animal.abd)=="Study"] <- "uniqueID"
## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
meta3 <- merge(animal.abd, ms.data, by="uniqueID")
meta3$last.grazing.event <- as.numeric(meta3$last.grazing.event)
meta3[,c("yi")] <- meta3[,c("yi")]
meta3[,"yr.grazed"] <- meta3$last.grazing.event
## compare climate & grazing duration variables against effect size
## get last year grazed as covariate
meta.ready <- subset(meta3, yr.grazed>0& yr.grazed < 60)
#mixed-effects meta-analysis for grazed vs ungrazed plots
m3 <- rma(yi=yi, vi=vi, mods=~ poly(yr.grazed,1), data = meta.ready, test="t")
summary(m3)
##
## Mixed-Effects Model (k = 107; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -113.7611 227.5222 233.5222 241.4841 233.7598
##
## tau^2 (estimated amount of residual heterogeneity): 0.2984 (SE = 0.0593)
## tau (square root of estimated tau^2 value): 0.5462
## I^2 (residual heterogeneity / unaccounted variability): 86.80%
## H^2 (unaccounted variability / sampling variability): 7.58
## R^2 (amount of heterogeneity accounted for): 9.95%
##
## Test for Residual Heterogeneity:
## QE(df = 105) = 620.3301, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 105) = 7.0875, p-val = 0.0090
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.1458 0.0652 2.2352 0.0275 0.0165 0.2751 *
## poly(yr.grazed, 1) 1.8068 0.6787 2.6622 0.0090 0.4611 3.1524 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## generate column for FG
meta.ready[,"trophic"] <- gsub("^.*?ertebrate","",meta.ready$FG)
## line
preds <- predict(m3, yr.grazed=meta.ready$yr.grazed)
# cbPalette <- c("#77564C", "#4381C1", "#9AD4D6", "#4FB286", "#564787", "#08415C")
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#D55E00", "#0072B2","#000000")
## animals
plot1 <- ggplot(meta.ready, aes(x=yr.grazed, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2)+ geom_jitter(aes(color=trophic, shape=higher.taxa), size=4, width=1) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Time since exclusion (yrs)") + ylab("Effect of livestock exclusion (LRR)") + geom_line(data=data.frame(preds), aes(x=meta.ready$yr.grazed, y=pred), lwd=1) +
geom_ribbon(data=data.frame(preds), aes(x=meta.ready$yr.grazed, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) + scale_color_manual(values=cbPalette)+
annotate("text", x=32, y=2.4, label="animal abundance", size=8)
## Warning: Ignoring unknown aesthetics: y
#### diversity
## Load in effect sizes
occurdat<-list.files( pattern="Eff.csv$",full=T)
animal.abd <- read.csv(occurdat[2])
colnames(animal.abd)[colnames(animal.abd)=="Study"] <- "uniqueID"
## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.revised.csv")
meta3 <- merge(animal.abd, ms.data, by="uniqueID")
meta3$last.grazing.event <- as.numeric(meta3$last.grazing.event)
meta3[,c("yi")] <- meta3[,c("yi")]
meta3[,"yr.grazed"] <- meta3$last.grazing.event
## compare climate & grazing duration variables against effect size
## get last year grazed as covariate
meta.ready <- subset(meta3, yr.grazed>0& yr.grazed < 60)
#mixed-effects meta-analysis for grazed vs ungrazed plots
m3 <- rma(yi=yi, vi=vi, mods=~ poly(yr.grazed,2), data = meta.ready, test="t")
summary(m3)
##
## Mixed-Effects Model (k = 26; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -1.7467 3.4934 11.4934 16.0353 13.7156
##
## tau^2 (estimated amount of residual heterogeneity): 0.0340 (SE = 0.0172)
## tau (square root of estimated tau^2 value): 0.1843
## I^2 (residual heterogeneity / unaccounted variability): 64.97%
## H^2 (unaccounted variability / sampling variability): 2.85
## R^2 (amount of heterogeneity accounted for): 13.85%
##
## Test for Residual Heterogeneity:
## QE(df = 23) = 64.9276, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 23) = 2.0041, p-val = 0.1576
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.1981 0.0492 4.0299 0.0005 0.0964 0.2999 ***
## poly(yr.grazed, 2)1 -0.0935 0.2567 -0.3642 0.7191 -0.6244 0.4375
## poly(yr.grazed, 2)2 -0.5362 0.2679 -2.0020 0.0572 -1.0903 0.0178 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## generate column for FG
meta.ready[,"trophic"] <- gsub("^.*?ertebrate","",meta.ready$FG)
## line
preds <- predict(m3, yr.grazed=meta.ready$yr.grazed)
# cbPalette <- c("#77564C", "#4381C1", "#9AD4D6", "#4FB286", "#564787", "#08415C")
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#D55E00", "#0072B2","#000000")
## animals
plot2 <- ggplot(meta.ready, aes(x=yr.grazed, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2)+ geom_jitter(aes(color=trophic, shape=higher.taxa), size=4, width=1) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Time since exclusion (yrs)") + ylab("Effect of livestock exclusion (LRR)") + geom_line(data=data.frame(preds), aes(x=meta.ready$yr.grazed, y=pred), lwd=1) +
geom_ribbon(data=data.frame(preds), aes(x=meta.ready$yr.grazed, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) + scale_color_manual(values=cbPalette)+
annotate("text", x=32, y=1, label="animal diversity", size=8)
## Warning: Ignoring unknown aesthetics: y
### Old figure
# ### calculate predicted risk ratios for full grazing time gradient
# preds <- predict(m3, newmods=c(0:35), transf=exp)
#
# ### radius of points will be proportional to the inverse standard errors
# ### hence the area of the points will be proportional to inverse variances
# size <- 1 / sqrt(meta.ready$vi)
# size <- size / max(size) *0.95
#
# ### set up plot (risk ratios on y-axis, absolute temperature on x-axis)
# plot(NA, NA, xlim=c(0,35), ylim=c(0.5,5),
# xlab="Year since control was last grazed", ylab="Effect of grazing relative to control (LRR)",
# las=1, bty="l", log="y", cex.lab=1.5, cex.axis=1.2)
#
# ### add points
# col <- data.frame(FG=unique(meta.ready$FG), colors= c("#999999","#E69F00", "#56B4E9", "#D55E00","#0072B2","#E69F00","#0072B2"))
# ptshape <- data.frame(FG=unique(meta.ready$FG), shape= c(21,21,21,21,21,23,23))
# meta.ready <- merge(meta.ready, col, by="FG")
# meta.ready <- merge(meta.ready, ptshape, by="FG")
# points(meta.ready$yr.grazed, exp(meta.ready$yi), pch=meta.ready$shape, bg=as.character(meta.ready$colors), cex=2)
# arrows(meta.ready$yr.grazed, exp(meta.ready$yi)-meta.ready$vi, ## lower
# meta.ready$yr.grazed, exp(meta.ready$yi)+meta.ready$vi, ## upper
# length=0.05, angle=90, code=3)
#
# # symbols(meta.ready$yr.grazed, exp(meta.ready$yi), circles=size, inches=FALSE, add=TRUE, bg=as.character(meta.ready$colors))
#
# legend(0,5.5, legend=col$FG, pt.bg=as.character(col$colors), pch=meta.ready$shape, cex=0.8)
#
# ### add predicted values (and corresponding CI bounds)
# lines(0:35, preds$pred, lwd=2)
# lines(0:35, preds$ci.lb, lty="dashed")
# lines(0:35, preds$ci.ub, lty="dashed")
### Load main data
meta <- read.csv("data//binary.simplified.csv")
## load studies used
div <- read.csv("diversityEff.csv")
abd <- read.csv("abundanceEff.csv")
studiesUsed <- c(div$Study, abd$Study)
dataused <- meta %>% filter(uniqueID %in% studiesUsed)
## count taxa
taxa <- dataused %>% group_by(Taxa) %>% summarize(n=length(unique(uniqueID)))
taxa
## # A tibble: 17 x 2
## Taxa n
## <fct> <int>
## 1 "" 4
## 2 "AMF" 1
## 3 "Amphibian" 2
## 4 "Amphibians" 2
## 5 "Annelida" 1
## 6 "Arthropoda" 41
## 7 "Aves" 12
## 8 "Insecta" 1
## 9 "Mammalia" 31
## 10 "Microbes" 1
## 11 "Mollusca" 1
## 12 "Nematoda" 2
## 13 "Nematode" 1
## 14 "Orthoptera" 1
## 15 "plants" 1
## 16 "Plants" 27
## 17 "Reptilia" 8
## unique studies
length(unique(studiesUsed))
## [1] 90
## list species
dataused %>% group_by(Detailed.Taxa ) %>% summarize(n=length(unique(Detailed.Taxa ))) %>% data.frame(.)
## Detailed.Taxa n
## 1 1
## 2 Acari 1
## 3 Amphibians 1
## 4 Annelida 1
## 5 Arachnida 1
## 6 Araneae 1
## 7 Arthropoda 1
## 8 Bee 1
## 9 Bovidae 1
## 10 Coleoptera 1
## 11 Collembola 1
## 12 Dasyuromorphia 1
## 13 Diptera 1
## 14 Equidae 1
## 15 Formicidae 1
## 16 Gastropoda 1
## 17 Hemiptera 1
## 18 Hymenoptera 1
## 19 Isoptera 1
## 20 Lagomorpha 1
## 21 Lepidoptera 1
## 22 Mycorrhiza 1
## 23 Orthoptera 1
## 24 plant 1
## 25 Psocoptera 1
## 26 Rodenta 1
## 27 Squamata 1
## 28 Thysanoptera 1
## load studies used
div <- read.csv("diversityEff.csv")
abd <- read.csv("abundanceEff.csv")
plt <- read.csv("plantEff.csv")
## merge data
plt <- subset(plt, measure=="diversity")
## select matching studies
abd <- abd[abd$Study %in% plt$Study,]
plt <- plt[plt$Study %in% abd$Study,]
## generate dataframe
plt <- data.frame(Study=plt$Study, plteff = plt$yi)
pltabd <- merge(plt, abd, by="Study", all.x=T)
## Model
m1 <- rma(yi=yi, vi=vi, mods=~ plteff, data = pltabd, test="t")
summary(m1)
##
## Mixed-Effects Model (k = 15; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -6.5583 13.1166 19.1166 20.8114 21.7832
##
## tau^2 (estimated amount of residual heterogeneity): 0 (SE = 0.0234)
## tau (square root of estimated tau^2 value): 0
## I^2 (residual heterogeneity / unaccounted variability): 0.00%
## H^2 (unaccounted variability / sampling variability): 1.00
## R^2 (amount of heterogeneity accounted for): 100.00%
##
## Test for Residual Heterogeneity:
## QE(df = 13) = 12.5828, p-val = 0.4805
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 13) = 13.2080, p-val = 0.0030
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.3604 0.0730 4.9378 0.0003 0.2027 0.5181 ***
## plteff 1.1351 0.3123 3.6343 0.0030 0.4603 1.8098 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## separate trophic group
pltabd[,"trophic"] <- gsub("^.*?ertebrate","",pltabd$FG)
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#D55E00", "#0072B2","#000000")
## predict line
preds <- predict(m1, newdata=data.frame(plteff=pltabd$plteff))
##plot
ggplot(data=pltabd, aes(x=plteff, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2)+geom_vline(xintercept=0, lwd=1, lty=2)+ geom_point(aes(color=trophic, shape=higher.taxa), size=4) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Effect of livestock exclusion on plant diversity (LRR)") + ylab("Effect of livestock exclusion on animal abundance (LRR)") + geom_line(data=data.frame(preds), aes(x=pltabd$plteff, y=pred), lwd=1) +
geom_ribbon(data=data.frame(preds), aes(x=pltabd$plteff, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) +
scale_color_manual(values=cbPalette)
## Warning: Ignoring unknown aesthetics: y
# annotate("text", x=-.25, y=-1.6, label="plant diversity \n animal abundance", size=6)
### Animal diversity
## merge data
## select matching studies
div <- div[div$Study %in% plt$Study,]
plt <- plt[plt$Study %in% div$Study,]
## generate dataframe
#plt <- data.frame(Study=plt$Study, plteff = plt$yi)
pltdiv <- merge(plt, div, by="Study", all.x=T)
## Model
m1 <- rma(yi=yi, vi=vi, mods=~ poly(plteff,1), data = pltdiv, test="t")
summary(m1)
##
## Mixed-Effects Model (k = 8; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 1.9345 -3.8691 2.1309 1.5062 14.1309
##
## tau^2 (estimated amount of residual heterogeneity): 0.0098 (SE = 0.0158)
## tau (square root of estimated tau^2 value): 0.0992
## I^2 (residual heterogeneity / unaccounted variability): 35.87%
## H^2 (unaccounted variability / sampling variability): 1.56
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 6) = 8.6282, p-val = 0.1956
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 6) = 0.0059, p-val = 0.9413
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## intrcpt 0.2353 0.0604 3.8940 0.0080 0.0874 0.3831 **
## poly(plteff, 1) -0.0117 0.1524 -0.0768 0.9413 -0.3847 0.3612
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## separate trophic group
pltdiv[,"trophic"] <- gsub("^.*?ertebrate","",pltdiv$FG)
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#D55E00", "#0072B2","#000000")
## predict line
preds <- predict(m1, plteff=pltdiv$plteff)
##plot
ggplot(data=pltdiv, aes(x=plteff, y=yi))+ geom_hline(yintercept=0, lwd=1, lty=2) + geom_vline(xintercept=0, lwd=1, lty=2)+ geom_point(aes(color=trophic, shape=higher.taxa), size=4) + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), text = element_text(size=16)) + xlab("Effect of livestock on plants (LRR)") + ylab("Effect of livestock on animals (LRR)") +# geom_line(data=data.frame(preds), aes(x=pltdiv$plteff, y=pred), lwd=1) +
#geom_ribbon(data=data.frame(preds), aes(x=pltdiv$plteff, y=pred, ymin=ci.lb, ymax=ci.ub), alpha=0.2) +
scale_color_manual(values=cbPalette)+ ylim(-0.5,0.25)+
annotate("text", x=-.3, y=-.4, label="plant diversity \n animal diversity", size=6)
## Warning: Removed 3 rows containing missing values (geom_point).